\chapter{时间序列分析}

\section{阈值非线性：时间序列}
加载\texttt{tsDyn}包，该包对普通线性模型\texttt{AR(m)}的表述为：
\begin{equation}\label{ts1}
 x_{t+s}=\phi_0+\phi_1x_{t}+\phi_2x_{t-d}+\phi_3x_{t-2d}+\cdots+\phi_mx_{t-(m-1)d} +\varepsilon_{t+s}
\end{equation}

实际上，$ d=1 $就是我们经常看到的\texttt{AR}模型形式。
\subsection{SETAR：自激励}
\[x_{t+s}=\begin{cases}
\phi_{1}+\phi_{10}x_t+\phi_{11}x_{t-d}+\cdots+\phi_{1L}x_{t-(L-1)d}+\varepsilon_{t+s}\quad Z_t\le th\\
\phi_{2}+\phi_{20}x_t+\phi_{21}x_{t-d}+\cdots+\phi_{2L}x_{t-(H-1)d}+\varepsilon_{t+s}\quad Z_t> th
\end{cases}\]

其中，$ Z_t $是阈值变量，$ th $是阈值。
\begin{itemize}
	\item 若阈值变量是本身的滞后值，即$ Z_t=x_{t-\delta d} $那么可以输入：\\
	\lstinline|obj <- setar(x, m=, d=, steps=, thDelay= )|
	
	其中，\lstinline|m,d|见\eqref{ts1}式， \lstinline|steps|是预测的步数$ s $，\lstinline|thDelay|是$\delta $。
	\item 若阈值变量是外生变量，则可以输入：\\
	\lstinline|obj <- setar(x, m=, d=, steps=, thVar= )|
	\item 若包含2个阈值，即3个体制，则可以输入：\\
	\lstinline|obj <- setar(x, m=, d=, steps=, thDelay=, nthresh=2)|
	\item 若想更改不同区制下的模型阶数设定，则：\\
	\lstinline|obj <- setar(x, d=, steps=, thDelay = , mL =, mH =)|
\end{itemize}

\paragraph{股指与通胀}我们有一个直觉：不同通胀水平下，股指的变化是不一样的。通胀较高时，股指似乎更容易涨。
\begin{figure}[H]
	\begin{lstlisting}
	rm(list = ls())
	setwd('E:\\17_Huadong\\teach\\MyLecture\\R')
	load('SH_CPI.RData')
	library(tsDyn)
	library(xts)
	library(magrittr)
	SH_CPI <- xts(SH_CPI[,-1],order.by = SH_CPI$DATETIME)
	# get yield
	SH_CPI$SH <- log(SH_CPI$M0020188.y) %>% diff()
	# lag CPI 
	SH_CPI$lagCPI <- lag(SH_CPI$M0000705.y)
	# delete 1st row that is NA
	DSH_CPI <- SH_CPI[-1,]
	# CPI have effect for SH in 2 month
	setar(DSH_CPI$SH,m=1,thVar = DSH_CPI$lagCPI) %>% summary()
	setar(DSH_CPI$SH,m=2,thVar = DSH_CPI$lagCPI) %>% summary()
	\end{lstlisting}
\end{figure}

二阶滞后的估计结果如下：
\begin{figure}[H]
	\includegraphics[scale=.8]{ThRlt.png}
\end{figure}

\subsection{LSTAR：平滑转移}
\begin{align*}
x_{t+s}=&(\phi_{1}+\phi_{10}x_t+\phi_{11}x_{t-d}+\cdots+\phi_{1L}x_{t-(L-1)d})(1-G(Z_t,\gamma,th))+\\
&(\phi_{2}+\phi_{20}x_t+\phi_{21}x_{t-d}+\cdots+\phi_{2L}x_{t-(H-1)d})G(Z_t,\gamma,th)+\varepsilon_{t+s}
\end{align*}

其中$ G $是Logit函数,
\[ G(Z_t,\gamma,th)=\frac{1}{1+e^{-\gamma(Z_t-th)}} \]

可以看到，$ Z_t\rightarrow -\infty, G\rightarrow 0; Z_t\rightarrow +\infty, G\rightarrow 1 $。即阈值变量的变化使得系数的改变是平滑的，而非陡峭的。图\ref{f1}的绘制代码在第\ref{chpDraw}章给出。
\begin{figure}[H]
	\includegraphics[scale=.6]{LSTAR}
	\caption{不同$ \gamma $下阈值变量$ Z $与$ G $的变化关系}\label{f1}
\end{figure}

估计的命令为\lstinline|lstar|，格式同\lstinline|setar|，只是要多设置一个\lstinline|gamma|。
\subsection{TVAR：阈值向量自回归}
考虑长期利率与短期利率之间的互动，构建一个2阶滞后的阈值向量自回归：

当$ long_{t-1}>\gamma $时，
\begin{align*}
long_t&=\beta_{10}+\beta_{11}^1long_{t-1}+\beta_{12}^1long_{t-2}+\beta_{11}^2short_{t-1}+\beta_{12}^2short_{t-2}+\varepsilon_{1t}\\
short_t&=\beta_{20}+\beta_{21}^1long_{t-1}+\beta_{22}^1long_{t-2}+\beta_{21}^2short_{t-1}+\beta_{22}^2short_{t-2}+\varepsilon_{2t}
\end{align*}

当$ long_{t-1}\le\gamma $时，
\begin{align*}
long_t&=\alpha_{10}+\alpha_{11}^1long_{t-1}+\alpha_{12}^1long_{t-2}+\alpha_{11}^2short_{t-1}+\alpha_{12}^2short_{t-2}+\varepsilon_t\\
short_t&=\alpha_{20}+\alpha_{21}^1long_{t-1}+\alpha_{22}^1long_{t-2}+\alpha_{21}^2short_{t-1}+\alpha_{22}^2short_{t-2}+\varepsilon_t
\end{align*}

估计代码为，
\begin{lstlisting}
library(tsDyn)
data(zeroyld)
TVAR(zeroyld, lag=2, nthresh=1, thDelay=1, trim=0.1, mTh=1, plot=TRUE)
\end{lstlisting}

其中，\texttt{lag}就是VAR的滞后，\texttt{nthresh}是阈值的数目，\texttt{thDelay}是阈值是滞后几阶的，\texttt{mTh}是哪个作阈值，第一个还是第二个？
\subsection{阈值非线性：面板数据}
考虑教育水平对城镇化的影响，可能较低的教育水平对城镇化影响较低，较高的教育水平对城镇化的影响较大。基于该思路，利用全国31个省的面板数据，估计一个面板阈值模型。模型设置如下，
\begin{align*}
UR_{ij}=&\beta_1+\beta_{14}HSC_{it}I(HSC_{it}\le \gamma)+\beta_{24}HSC_{it}I(HSC_{it}> \gamma)+\beta_{5}IGDP_{it}+\\
& \beta_{6}URIG_{it}+\beta_{7}EO_{it}+\beta_{8}IP_{it}+\varepsilon_{it}
\end{align*}

其中，$UR$是城镇化水平，$ HSC $是人力资源水平。其他是控制变量。估计面板阈值，建议调用stata程序，R语言中的\texttt{pdR}包不是很好。

\begin{figure}[H]
	\begin{lstlisting}
	load('data180201.RData')
	names(data01900) <- c('region','year','CityRatio','base','med','high','gdp','dist',
	'ratio','open','cityRatio1','human','EduFee')
	data01900$region <- as.factor(data01900$region)
	library(RStata)
	
	StataCode <- '
	set more off
	xtset region year
	xthreg CityRatio gdp dist ratio open, rx(human) qx(human) thnum(1) trim(0.01) bs(300)
	'
	stata(StataCode, stata.version = 14,data.in = data01900)
	\end{lstlisting}
\end{figure}

其中，\texttt{rx}是体制依赖自变量，\texttt{qx}是阈值变量，\texttt{thnum}是阈值个数，\texttt{bs}是自助抽样次数。

\section{Kalman滤波: \texttt{MARSS}}
\subsection{理论模型}
当存在某一个变量，它是随时间而变化的，一般而言，该变量的真值我们无法得到，但我们往往能得到该变量经过某个黑箱之后的观测值。具备这种特征的变量(即后文所述的状态变量)，我们就能够使用卡尔曼滤波系统为其建模。

经典的理论形式：
\begin{align}\label{klm1}
\bm{x}_t&=\bm{u}_t+\bm{B}_t\bm{x}_{t-1}+\bm{w}_t,\; \bm{w}_t\sim MVN(0,\bm{Q}_t)\\\label{klm2}
\bm{y}_t&=\bm{a}_t+\bm{Z}_t\bm{x}_{t}+\bm{v}_t,\; \bm{v}_t\sim MVN(0,\bm{R}_t)\\\nonumber
\bm{x}_0& \sim MVN(\bm{\pi},\bm{V})
\end{align}

其中，$ \bm{y}_t $就是可观测变量，$ \bm{x}_t $就是状态变量。

通俗地说，卡尔曼滤波本质上是对一个观测不到的状态变量$x_t$进行建模，虽然我们观测不到该状态变量，但我们一方面知道该状态变量的演进形式，即状态方程\eqref{klm1}式，一方面我们还知道该状态变量在经过某个黑箱（该黑箱实际上就是度量方程\eqref{klm2}式）后，我们能够得到该状态变量的一个观测值，即$ \bm{y}_t $。于是，通过一定的估计方法，我们可以获得状态方程和度量方程中的参数估计，这样在状态方程和度量方程已知的情况下，通过输入状态变量的初值$ \bm{x}_0 $，我们就获得了所有的时期观测不到的状态变量的取值。总之，我们所做的这一切目的就是要把观测不到的状态变量通过度量方程和状态方程将它算出来，显化出来。




\subsection{Kalman滤波：一个例子}
注意数据的排列，对于观测变量$ \bm{y} $，其维数为$ n\times T $, 行为观测的序列个数，列为时间步数。

在估计中，要设定的参数就包括$B,u,Q,Z,A,R,x0,V0 $。设定参数的一个原则就是，若是字母则需要估计，若是数字则不需要估计。

